source('codes/00_import_libraries.R')


#-----------------------------------------------------------------------------#
# Data ----
#-----------------------------------------------------------------------------#

dict_name = 'base'
macro_list = c('DP', 'EP', 'SVAR', 'TBL')
gw = readRDS('input/gw_macro_predictors.RDS')$dt
dt = readRDS('input/combined_data.rds')[[dict_name]] %>% 
  select(-OG) %>% 
  mutate(
    ret  = log(1 + ret0),
    ret1 = log(1 + ret1),
    e_War = arima(War, order = c(1, 0, 0)) %>% resid()
  ) %>% 
  left_join(gw, by = 'ym') %>% 
  select(ym, ret, ret1, War, e_War, PLS_All, PLS_All_recursive, any_of(macro_list)) %>% 
  drop_na(ret)



#-----------------------------------------------------------------------------#
# VAR decomposition function ----
#-----------------------------------------------------------------------------#

fn_CF_DR = function(macro_list, xvar = 'War') {
  
  #---------------------------------------#
  # VAR decomposition ----
  #---------------------------------------#
  
  y_list = c('ret', macro_list)
  
  n = nrow(dt)
  k = length(y_list)
  Y1 = dt[2:n, y_list] %>% as.matrix()
  Y0 = dt[1:(n-1), y_list] %>% as.matrix()
  
  # VAR: Y1 = A.Y0 + U1
  # coefficient matrix: each column for each equation
  A = solve(t(Y0) %*% Y0) %*% (t(Y0) %*% Y1)
  
  # residuals
  U1 = Y1 - Y0 %*% A
  
  # return innovation
  e1 = c(1, rep(0, k-1))
  e_R = U1 %*% e1
  
  # DR and CF news
  rho = 0.96
  e_DR = U1 %*% (A %*% solve(diag(k) - rho*A)) %*% e1
  e_CF = e_R + e_DR
  
  # expected returns
  R_hat = Y0 %*% A %*% e1
  
  
  #---------------------------------------#
  # Predict with War ----
  #---------------------------------------#
  
  # combine with War
  df = cbind(dt[1:(n-1), c('ym', xvar, 'ret1')], R_hat, e_CF, e_DR) %>% 
    as_tibble() %>% 
    drop_na() %>% 
    mutate(
      across(xvar, ~ (.x - mean(.x))/sd(.x)),
      # War = (War - mean(War))/sd(War),
      across(c(ret1, R_hat, e_CF, e_DR), ~.x*12*100)
    )
  
  
  yvars = c('ret1', 'R_hat', 'e_CF', 'e_DR')
  lapply(yvars, function(y) {
    fml = glue('{y} ~ {xvar}') %>% as.formula()
    mod = lm(fml, df)
    mod_robust = coeftest(mod, vcov. = NeweyWest(mod, lag = 1, prewhite = 0))
    
    out = huxreg(
      mod_robust,
      coefs           = xvar,
      error_format    = '({statistic})',
      stars           = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
      statistics      = c(0),
      borders         = 0,
      outer_borders   = 0,
      note            = NULL
    )
    
    return(out[-1, -1])
  }) %>% 
    do.call(what = cbind)
}


#-----------------------------------------------------------------------------#
# Apply ----
#-----------------------------------------------------------------------------#

macro_list = list(
  'DP',
  c('DP', 'EP'),
  c('DP', 'SVAR'),
  c('DP', 'TBL'),
  c('DP', 'EP', 'SVAR', 'TBL')
)

macro_names = sapply(macro_list, function(x) paste0(c('R', x), collapse = ', '))
row_names = lapply(macro_names, function(x) paste0(c(x, ''))) %>% unlist()
col_names = c('', '$r_{t+1}$', '$E_t r_{t+1}$', '$CF_{t+1}$', '$DR_{t+1}$')

xvar = 'War'
output = lapply(macro_list, \(z) {
  fn_CF_DR(macro_list = z, xvar = xvar)
}) %>% 
  do.call(what = rbind) %>% 
  set_number_format(2) %>% 
  insert_column(row_names) %>% 
  insert_row(col_names) %>% 
  set_align(1, everywhere, 'center') %>% 
  set_tb_borders(1, everywhere, 0.8) %>% 
  set_bottom_border(final(1), everywhere, 0.8) %>% 
  set_all_padding(0) %>% 
  set_escape_contents(FALSE)



#------------------------------------------------------------------------
# Export results
#------------------------------------------------------------------------

quick_xlsx(output, file = glue('output/tables/Table_11.xlsx'))
# cat(to_latex(output, tabular_only = TRUE), file = glue('output/tables/Table_11.xlsx'))












